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Abstract 

We study the thermodynamic Casimir force for films in the three-dimensional Ising universality 
class with symmetry breaking boundary conditions. We focus on the effect of corrections to scaling 
and probe numerically the universality of our results. In particular we check our hypothesis that 
corrections are well described by an effective thickness Lo,e// = Lo+c(Lo+L s ) 1- ^ +L s , where c and 
L s are system specific parameters and u ~ 0.8 is the exponent of the leading bulk correction. We 
simulate the improved Blume-Capel model and the spin-1/2 Ising model on the simple cubic lattice. 
Note that in contrast to the Ising model, for the improved Blume-Capel model the amplitude of 
corrections oc Lq w is strongly suppressed. First we analyse the behaviour of various quantities 
at the critical point. Taking into account corrections oc Lq w in the case of the Ising model, 
we find good consistency of results obtained from these two different models. In particular we 
get from the analysis of our data for the Ising model for the difference of Casimir amplitudes 

A H — A_|__|_ = 3.200(5), which nicely compares with A^ — A + _|_ = 3.208(5) obtained by studying 

the improved Blume-Capel model. Next we study the behaviour of the thermodynamic Casimir 
force for large values of the scaling variable x = i[Lo/£o]- This behaviour can be obtained up to 
an overall amplitude by expressing the partition function of the film in terms of eigenvalues and 
eigenstates of the transfermatrix and boundary states. Here we show how this overall amplitude 
can be computed with high accuracy. We find good agreement of the numerical results obtained by 
studying the spin-1/2 Ising and the improved Blume-Capel model. Finally we discuss our results 

for the scaling functions 8^ and 9 ++ of the thermodynamic Casimir force for the whole range of 

the scaling variable. We conclude that our numerical results are in accordance with universality. 
Corrections to scaling are well approximated by an effective thickness. 

PACS numbers: 05.50.+q, 05.70.Jk, 05.10.Ln, 68.15. +c 
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I. INTRODUCTION 



At a second order phase transition various quantities like the correlation length £ or the 
specific heat Cbuik diverge following power laws such as 

£~£ ,±|*r' C bulk ~A ± \t\- a (1) 

where t — (T — T c )/T c is the reduced temperature, v and a are the critical exponents of the 
correlation length and the specific heat, respectively. The indices ± of the amplitudes £ ,± 
and A± indicate the phase: + for the high temperature phase and — for the low temperature 
phase. Critical exponents such as v and a and amplitude ratios such as £ ,+/£o,- an d A + /A_ 
are universal. This means that these quantities do not depend on the microscopic details of 
the system but are exactly the same for all systems within a universality class. A universality 
class is characterized by the dimension of the system, the range of the interaction and the 
symmetry properties of the order parameter. Note that it is a non-trivial task to identify 
an order parameter. In particular, the symmetry properties of the order parameter are not 
necessarily the same as those of the Hamiltonian. For reviews on critical phenomena see e.g. 

Power laws such as eq. (CQ) are valid only asymptotically in the limit t — > 0. At finite 
reduced temperature corrections have to be taken into account [H @] 



£ = f ,±|*r" x (l + a±\t\ e + bt + c±\t\ 26 + d±\tf + . 



(2) 



There are analytic and non-analytic (confluent) corrections. The non-analytic corrections are 
associated with non-trivial exponents 9 = uu, 6' = vu' , ... . For the universality class of the 
three-dimensional Ising model one finds consistently u ~ 0.8 from field theoretic methods, 
the analysis of high temperature series expansions and Monte Carlo simulations of lattice 
models Q. Our recent estimate is u> = 0.832(6) [7]. The estimate u' = 1.67(11) obtained 
by the scaling field method [s| still lacks confirmation by other approaches. Furthermore we 
expect corrections caused by the breaking of symmetries by the lattice. In the case of the 
simple cubic lattice that we consider here, these corrections are associated with to" ~ 2 [§]. 

The singular behaviour ([I]) requires that the thermodynamic limit is taken. For finite 
systems, the behaviour of thermodynamic quantities is given by analytic functions of the 
parameters of the system and its linear size Lq. Finite size scaling [Io[ predicts that in the 
neighbourhood of the critical point, for sufficiently large Lq, this behaviour is characterized 
by a universal function of certain combinations of the parameters of the system and its 
linear size Lq. In the absence of an external field, a quantity A(L Q ,t) that is a function of 
the temperature and the linear size L Q of the system behaves as 

A(L ,t)^L^(t[L /6 )+ ] 1/i/ ) (3) 

where the function g(x) depends on the universality class of the bulk system and on the 
geometry of the finite system and y = w/u, where A(oo,t) oc \t\~ w . Also finite size scaling 



is affected by corrections to scaling [10 



A(L , t) = LI s(t[L o /U+] lA [1 + b q{t[LQ/iQ^) LJ* + ...] (4) 
where q(x) is a universal function and b depends on the details of the system. 
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Here we shall study films with symmetry-breaking boundary conditions. In addition to 
the corrections discussed above, these boundary conditions give rise to additional corrections, 
where the leading one is oc Lq 

pea, 

where Lq is now the thickness of the film. In this 
work we focus on the thermodynamic Casimir effect [l4j in films. Due to the fact that in 
the neighbourhood of the critical point the range of thermal fluctuations is restricted by the 
finite thickness of the film an effective force arises. The thermodynamic Casimir force per 
area is given by 

r Casimir o j j yJ) 

where f ex = ffu m — L Q f bu i k is the excess free energy per area of the film, where ffu m is the 
free energy per area of the film and fbuik the free energy density of the bulk system. The 
thermodynamic Casimir force per area follows the finite size scaling law 



Casimir 



k B TLf 9{t[L Q /^ + ]^) 



(6) 



see e.g. ref. [15[. After the seminal work [14] it took about two decades until the thermo- 
dynamic Casimir effect could be demonstrated in experiments. The data obtained for films 
of different thicknesses of 4 He near the A-transition are represented to a reasonable approxi- 
mation by a unique finite size scaling function [l6|, 13] . Also experiments with liquid binary 
mixtures near the mixing-demixing were performed, where either films 18] or the sphere- 
plate geometry [19|, |20| were studied. Unfortunately, field theoretic methods do not allow to 



compute the scaling function 9{x) for the full range of the scaling variable [12J, |13J. There 



fore it was an important achievement that recently the thermodynamic Casimir force was 
computed by Monte Carlo simulations of lattice models. Correspondin g to the experiments 
on 4 He, the XY model on the simple cubic lattice was simulated [H], H[. Also the Ising 
model on the simple cubic lattice that shares the universality class of the mixing-demixing 
transition of binary mixtures was studied 22|, |23|. A reasonable match of the universal scal- 
ing functions 9 obtained from experiments and the corresponding Monte Carlo simulations 
of lattice models was found. For a recent review see [lij . 

However it turned out that it is quite difficult to obtain precise results for the universal 
scaling function 9 from these Monte Carlo simulations. For the thicknesses that can be 
reached, corrections to scaling are still significant. Fitting the data it is difficult to disen- 
tangle corrections oc and oc L^ 1 . Furthermore the universal function q{x), eq. (jlj), that 



governs the corrections oc L u is a priori unknown. The authors of |2lH23l] used ad hoc ap 



proximations of q(x) in the analysis of their data. Depending on the particular ansatz that 
they used, the results of [HI, 23[ for the universal scaling function vary by a large amount. 

In order to alleviate this problem we [25|, 26] studied improved models which are char- 
acterized by the fact that the amplitude of the leading bulk correction vanishes. Since the 
parameter of the improved model is determined numerically, in practice a residual amplitude 
remains, which is however at least by a factor of 30 smaller than that of the Ising model and 
the XY model on the simple cubic lattice, respectively [3, 13]. O ur results for the scaling 



functions of the thermodynamic Casimir force agree qualitatively with those of refs. [21M23 



However the numerical discrepancies are considerably larger than the errors that are quoted. 
In particular, the results obtained very recently in [28j from simulations of the Ising model 
by using the prefered ansatz of the authors, eqs. (17,18) of [28J, deviate clearly from those 
of [H|; See fig. 6 a of [28|; and from that of ]29i]; See fig. 6 b of_0]. For a discussion of 
this fact by the authors of 28], see the text on page 041605-9 of [28] starting about 20 lines 
below table II. 
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The aim of the present work is to reach a better understanding of corrections to scaling. 
This means that we intend to determine the function q(x) of eq. (jlj) for the thermodynamic 
Casimir force. Note that due to universality of the function q(x) our results might also be 
useful in the analysis of data obtained in experiments. Also here we start with an ansatz for 
q(x) which is motivated as follows. The corrections oc Lq 1 caused by the boundaries can be 
expressed by a constant shift in the thickness of the film. In equations such as eq. the 
thickness L is replaced by 

Lo, e ff = Lq + L s , (7) 

where L s depends on the details of the system but not on the observable. Here we shall 
probe the hypothesis that in an analogue way corrections oc Lq U can be taken into account 
by 

L , eff = L + c(L + L,) 1 - u + L a . (8) 

While renormalization group arguments suggest that eq. (J7J) is indeed exact, the general- 
ization is at best a good approximation. It is motivated by the fact that for the strongly 
symmetry breaking boundary conditions studied here fluctuations are suppressed in the 
neighbourhood of the boundaries. Hence the effect of corrections to scaling should be the 
largest close to the boundaries. Plugging eq. flS} into eq. (j3J), ignoring the correction oc Lq 1 
due to the boundary, we get 



L l g{x) x 1 + c 



A(L ,t) = (L + cLl^y g(t[(L + cL ^) /£,,+] ^) 

2ur 



xg'{x) 

y + - 



V" + o( V) ( 9 ) 



v g(x) 

where x = £[L /£o,+] ■ Hence our hypothesis (|Sj) results in 

q(x)=y + - 9 -^. (10) 
v g[x) 

The outline of the paper is the following: In section |TT] we define the models that we 
simulated and the observables that we measured. In section [TV] we study various quantities 
exactly at the critical point. Next, in section |V] we study the behaviour of the thermody- 
namic Casimir force for large values of the scaling variable x. To this end, we analyse the 
magnetisation profile near the boundary of the film and the correlation function of the bulk 
system. In section [VTJ we discuss our results for the scaling functions 9 ++ and 9 + _ in the 
full range of the scaling argument. Then we summarize and discuss our results. Finally in 
the appendix we discuss various results obtained for the bulk of the spin-1/2 Ising model. 



II. MODEL 

We study the Blume-Capel model on the simple cubic lattice. It is defined by the reduced 
Hamiltonian 

H=-0 Yl s ^y + D J2 s " > ( n ) 

<xy> x 

where the spin might assume the values s x 6 {—1,0,1}. x = (xq,xi,X2) denotes a site 
on the simple cubic lattice, where Xj G {1,2, and < xy > denotes a pair of nearest 

neighbours on the lattice. The inverse temperature is denoted by (3 = l/k B T. The partition 
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function is given by Z = Y2{ s } ex P( — where the sum runs over all spin configurations. 
The parameter D controls the density of vacancies s x = 0. In the limit D — > — oo vacancies 
are completely suppressed and hence the spin-1/2 Ising model is recovered. 

In d > 2 dimensions the model undergoes a continuous phase transition for — oo < 
D < D tr i at a C that depends on D. For D > D tTi the model undergoes a first order 
phase transition. The authors of [30| give for the three-dimensional simple cubic lattice 
Dtn = 2.0313(4). 

Numerically, using Monte Carlo simulations it has been shown that there is a point 
(D*,/3 C (D*)) on the line of second order phase transitions, where the amplitude of leading 
corrections to scaling vanishes. Our recent estimate is D* = 0.656(20) j7j. In 0] we simulated 
the model at D = 0.655 close to /3 C on lattices of a linear size up to L = 360. From a standard 
finite size scaling analysis of phenomenological couplings like the Binder cumulant we find 
/3 C (0.655) = 0.387721735(25). Furthermore the amplitude of leading corrections to scaling 
is at least by a factor of 30 smaller than for the spin-1/2 Ising model. As discussed in the 
appendix I A II we shall use /3 C = 0.22165462(2) as estimate of the inverse critical temperature 
of the spin-1/2 Ising model in the following. 



In [3jj we simulated the Blume-Capel model at D = 0.655 in the high temperature phase 
on lattices of the size L 3 with periodic boundary conditions in all directions and L ^ 10£ 
for 201 values of /3. For a few values of /3 we performed new simulations that reduced the 
statistical error considerably. In particular for /3 = 0.3872, which was our value closest to 
(3 C , we get £ 2 nd(0.3872) = 26.7013(15) for second moment correlation length now. Taking 
into account these new data we arrive at the slightly revised result 



The analogue result for the spin-1/2 Ising model is given in eq. ( lAlOj) in Appendix I A 21 

In the high temperature phase there is little difference between ^2nd and the exponential 
correlation length £ e ™_which is defined by the asymptotic decay of the two-point correlation 
function. Following [32|: 



for the thermodynamic limit of the three-dimensional system. Note that in the following £o 
always refers to £2^,0,+ • 

A. Film geometry and boundary conditions 

In the present work we study the thermodynamic Casimir effect for systems with film 
geometry. In the ideal case this means that the system has a finite thickness L , while in 
the other two directions the thermodynamic limit L\ , L2 — > 00 is taken. In our Monte Carlo 
simulations we shall study lattices with L$ -^1,-^2 and periodic boundary conditions in 
the 1 and 2 directions. Throughout we shall simulated lattices with L\ = L2 = L. 

In the direction we take symmetry breaking boundary conditions. A strong breaking 
of the symmetry is achieved by fixing the spins at the boundary to either —1 or 1. Here we 
shall put these fixed spins on the layers at xq = and at xq — Lq + 1. This means that Lq 
gives the number of layers with fluctuating spins. In the following we shall consider the two 
choices: 




6nd,o,+ = 0.2283(1) - 1.8 X {y - 0.63002) + 275 x (f3 c - 0.387721735) 
using t = (3 C — (3 as definition of the reduced temperature. 



(12) 



1.000200(3) 



(13) 
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• ++ boundary conditions: s x — 1 for all x with xq — or x — Lq + 1. 

• H — boundary conditions: s x = 1 for all x with Xo = and s x = — 1 for all x with 
x = L + I. 



B. Free energy, energy and specific heat 

For bulk systems we define the reduced free energy density as 



fi 



bulk 



LqL\L 2 



(14) 



This means that compared with the free energy density fbulk, a factor k^T is skipped. 
Correspondingly we define the energy density as the derivative of minus the reduced free 
energy density with respect to (3 



E, 



bulk 



1 d\nZ 

L Q LiL 2 d(3 L, 



and the specific heat 



a 



dE, 



bulk 



bulk 



d/3 L LiL 2 



E 



S X Sy 



E 



S X Sy 



\<x,y> / I \<x,y> 

In the case of films we consider the reduced free energy per area 



/ 



1 



InZ 



and the energy per area 



E = 



1 d\nZ 
L X L 2 d/3 L X L 



-MY. 



i<x,y> 



(15) 



(16) 



(17) 



(18) 



C. The magnetization profile of films 

The film is invariant under translations in the 1 and 2 direction of the lattice. Therefore 
the magnetization only depends on x and we can average over x-y and x 2 : 



m(x ) 



L 2 



(19) 



Since the film is symmetric for ++ boundary conditions and anti-symmetric for H — bound- 
ary conditions under reflections at the middle of the film, m(x ) = m(L — x + 1) for ++ 
boundary conditions and m(x ) = —m(L — x + 1) for H — boundary conditions. 
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D. The correlation length 



The exponential correlation length £ of the bulk system is defined by the decay of the 
slice-slice correlation function 

G(r) ~ cexp(-r/£) (20) 
for large distances r. The slice-slice correlation function is defined as 

G(r) = (S(x )S(x + r)) (21) 

where 

S'(xo) = t== V (s X0 , Xl , X2 - (m» (22) 

where (m) is the bulk magnetisation that vanishes in the high temperature phase, for a 
vanishing external field. 

For a detailed discussion of the second moment correlation length defined for films see 
section III C of 1261. 



E. Monte Carlo algorithms 

In the case of the Ising model we simulated the films with Lq < 68 using a local Metropolis 
algorithm and a multispin coding implementation. We used the same program, up to small 
modifications to implement the boundary conditions, as discussed in section 3 of ref . [33] . On 
one core of an Intel (tm) Xeon(tm) E5520 CPU running at 2.27 GHz the program achieves 
1.9 x 10 9 spin updates per second. This is about 100 times faster than on the fastest 
workstation that was available to us in 1993. Most simulations were performed on Quad- 
Core AMD Opteron(tm) 2378 CPUs running at 2.4 GHz. Here the program achieves 1.4 x 10 9 
spin updates per second on one core. In relation with section |V] we simulated films with ++ 
boundary conditions with Lq > 68. These were simulated by using a special version of the 
cluster algorithm as discussed ref. [26]. In the case of the Blume-Capel model we simulated 
the films using the same algorithms as discussed in section V of ref. [26[ . 

Mostly we simulated lattices with periodic boundary conditions in all directions with the 



single-cluster algorithm [34| in the case of the Ising model and a hybrid [35| of the local 
heat-bath and the single-cluster algorithm in the case of the Blume-Capel model. 

In all our simulations we used the Mersenne twister algorithm [36j as pseudo-random 
number generator. In total our simulations took about the equivalent of 50 years of CPU 
time on a single core of a Quad-Core AMD Opteron(tm) Processor 2378 CPUs running at 
2.4 GHz. 



III. FINITE SIZE SCALING AND CORRECTIONS TO SCALING 

The reduced excess free energy per area of a film is given by 

f ex (Lo, 0) = f(L , 0) ~ Lof Mk (P) . (23) 

In the reduced excess free energy the analytic bulk contribution cancels. Therefore it can 
be written as 

f ex (L ,P) = f ex , s (L ,P) + 2f r (P) (24) 
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where f eXjS is the singular part and f r is an analytic contribution due to the boundaries. 
In the absence of an external field, this contribution is the same for a boundary where all 
spins are fixed to +1 and one where all spins are fixed to —1. The free energy of a system is 
conserved under renormalization group transformations. Therefore the singular part of the 
reduced excess free energy behaves as 

f ex , s (Lo,(3) = L 2 H(t[L /^ }y,bL y \...) (25) 

where y% = 1/V and y\ = —u are the thermal and the leading irrelevant renormalization 
group exponent, respectively. Expanding the universal scaling function H(x,u,...) in u 
around u = we arrive at 



f ex<s (L , 13) = Lfh{x) [I + 6 p(x)L^ + ...] 



(26) 



where x = t[Lo/£o] 1//!/ and the leading correction is characterized by the universal function 
p(x). Taking minus the derivative with respect to Lq we get the thermodynamic Casimir 
force 



k h T 



Fc 'asimir — Lq 6{x) 



1 + 61 p(x 



h(x) 
6{x) L 



X 

UP(X) - -P'(X) ) Lq u + 



where 



x 



d{x) = 2h{x) - -h'{x) . 



(27) 



(28) 



Note that at the critical point 8(0) = 2h(0). In the literature h(0) is called Casimir amplitude 
and is denoted by A. Also note that 



e\o) 



i 

2- - 

v 



h'(0) 



(29) 



Taking minus the derivative with respect to we get 

h(x)p'(x) 



E ex (LQ,(3) = LQ 2 [LQ/Co} 1/u h'l 



x 



1 + 6 p(x) + 



h'(x) 



T-U) 



2/;(/3) • (30) 



IV. FINITE SIZE SCALING AT THE CRITICAL POINT 

First we study finite size scaling at the critical point, i.e. x = 

t^o/^o] 1 ^ = 0. To this 

end we analyse data for the free energy difference between films with H — and ++ boundary 
conditions, the energy density and the magnetisation profile for both types of boundary 
conditions. Finally we also consider the second moment correlation length for H — boundary 
conditions. 

For a given quantity at a given value of a; it is a trivial recast to express corrections to 
scaling in the form ([S]) . The non-trivial question that we investigate here is whether leading 
corrections in different quantities can be expressed by the same or at least similar effective 
thicknesses L e ff. 

In the ansaetze below we shall use in addition to eq. ([8]) 

L , e // = L + L s + c(L + L a )- U + d(L + L s )~ e (31) 
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in order to probe for the effect of subleading corrections. As discussed in the intoduction, 
there are infinitely many subleading corrections starting with e = 2ui m u>', 1+oj and u" ~ 2. 
Given the accuracy of our data, it is only possible to put one subleading correction in the 
ansatz. In the following we shall take either e = 1.664 or e = 2. Fitting with ansaetze that 
only approximate the behaviour of the data one has to be aware of systematical errors. In 
the literature it is often implicitly assumed that an acceptable xVd.o.f. means that such 
systematical errors are small and of a similar size or even smaller than the statistical errors 
of the fit parameters. However this is definitely not the case. The severity of the problem 
depends of course on the type of the approximation and the range of the data that are 
available. Below we shall see that the differences between results of fits with eq. OH]) and 
ones with eq. fl3T|) are e.g. five times larger than the statistical error. The error that we 
quote for final results is chosen such that both the results of fits with eq. OH]) and eq. Q3~T]) 
are covered. 



A. The difference of free energies per area between H — and ++ boundary condi- 
tions 



First we studied the difference 

D /)+ _ )++ (L 0) P) = /+_(L , P) - U+(L , P) , (32) 

where /+_ and / ++ are the reduced free energies for H — and ++ boundary conditions, 
respectively. In this difference the surface and the bulk contributions exactly cancel and 
therefore at the critical point 

D f!+ _ >++ (L , C ) ~ (A + _ - A ++ ) Lf , (33) 

where A + _ and A ++ are the Casimir amplitudes for H — and ++ boundary conditions, 



respectively. Similar to the case of periodic and anti-periodic boundary conditions [37], |38 
the ratio Z_ + /Z ++ of partition functions can be directly computed by using the cluster 
algorithm. To this end one determines for ++ boundary conditions the fraction of cluster 
decompositions where the two boundaries do not belong to the same cluster. These cluster 
decompositions would allow to update to H — boundary conditions. Since for H — boundary 
conditions the update to ++ boundary conditions is always allowed, the fraction discussed 
above is an estimate of Z_ + /Z ++ . 

Unfortunately, at the critical point, for L ^> Lq, the ratio Z \-/Z ++ is far too small to 

allow for an efficient sampling. Therefore we simulated in the high temperature phase at 
P = 0o such that Lo/£(/?o) ~ 6, where £ is the bulk correlation length. Here, for L = 4L , 

which we used in our simulations, the value of Z \-/Z ++ is a few percent. In order to 

get / + _ — / ++ at larger values of p, in particular at the critical point, we performed an 
integration of energy differences: 

D f ^ ++ (L , P) = D Ji+ _, ++ (L , p ) - / dp D E ^ ++ (L , P) , (34) 

where De,-\ ,++ = E+- —E ++ . We performed this integration numerically, using the trape- 
zoidal rule. To this end, we used at least 36 values of P between p Q and p c as nodes. For a 
detailed discussion of the corresponding Monte Carlo simulations see section |VT] below. In 



9 



TABLE I. We give the difference Df _| _| of the reduced free energies per area between H — and 

++ boundary conditions at our estimates of the inverse critical temperature, i.e. (3 = 0.22165462 
for the Ising model and f3 = 0.387721735 for the Blume-Capel model at D = 0.655. 



L 


Model 


D f,+-,+- 


14 


1 


0.01069953(37) 


15 


1 


0.00953606(25) 


16 


1 


0.00855417(15) 


17 


1 


0.00771682(12) 


24 


1 


0.00423239(15) 


32 


1 


0.002522796(50) 


34 




0.002258418(55) 


48 




0.00119288(10) 


64 




0.000693495(64) 


68 




0.000617863(63) 


16 


BC 


0.00999910(67) 


17 


BC 


0.00897065(65) 


32 


BC 


0.00279016(11) 


34 


BC 


0.00248788(11) 


68 


BC 


0.00065641(11) 



most cases we used the same data as discussed in section |VH Only for the Ising model at 
the thicknesses Lq = 24 and 48 and the Blume-Capel model at the thickness Lq = 68 we 
performed additional simulations. For an analytic integrand, the estimate obtained by using 
the trapezoidal rule behaves as 1(h) = 1(0) + ah 2 + 0(/i 4 ), where 1(0) is the integral to be 
computed and h is the step-size. We estimated the systematic error by computing I(2h), 
i.e. performing the integration with half of the available data points. The systematic 
error is then estimated by e = (I(2h) — 1(h)) /3. It turned out that the systematic error e is 
considerably larger than the rather small statistical error. Therefore, we extrapolated our 
result as 1(0) = 1(h) — [I(2h) — I(h)]/3+ 0(h A ). In the case of the Blume-Capel model and 
Lq = 34, where we simulated at 116 values of (3 between /3 and (3 C we checked the efficiency 
of the extrapolation by computing 1(h), I(2h) and I (Ah). We found agreement between 
1(h) - (I(2h) - 1(h))/ 3 and I(2h) - (I(Ah) - I(2h))/3 within the statistical error. In table 
Uwe summarized our numerical results for the critical point. 

We fitted the data obtained for the Ising model with the ansaetze 

D fi+ ^ ++ = A[L + L S + c(L + Ls) 1 -"]- 2 (35) 

and 

= A [Lq + L s + c(L + L s y-" + d(L + L s ) l ~T 2 , (36) 

where we set either e = 1.664 or e = 2. 

Fitting with the ansatz fl35|) . setting u = 0.832 we get for L 0tmin = 16 the result A = 
3.1987(9), c = 1.429(12), L s = 1.043(12) and x 2 /d-o.f.= 1.20. Note that all data with L > 
Lo,min are taken into account in the fit. Instead, taking oj = 0.826 we get A = 3.1995(9), 
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c = 1.367(11), L s = 1.100(15) and x 2 /d.o.f.= l - 2Q - This means that the estimate of A 
depends little on the value of co, while c and L s are quite sensitive to it. We redid these fits 
for evaluated at (3 = 0.2216546. The results change only by little. 

Next we fitted all data, i.e. L ,mm = 14, with the ansatz (1561) . We get, fixing oj = 0.832 
and e = 1.664 the results A ='3.2025(21), c = 1.57(6), L s = 0.78(11), d = 0.35(14) 
and x 2 /d-o.f.= 1.47. Instead, for e = 2 we get A = 3.2016(19), c = 1.52(4), L s = 0.88(7), 
d = 1.52(4) and x 2 /d.o.f.= 1.45. We see that by adding a subleading correction the estimate 
of A changes little, while the results for c and L s are considerably shifted. Note that 
the estimates of c and L s are highly anti-correlated. The resulting £o,e//, eq. ([8]), for the 
thicknesses analysed here, depend much less on the ansatz that is used. Taking all fits 
discussed above into account we conclude 

A + _ - A ++ = 3.200(5) . (37) 

Next we fitted our data for the Blume-Capel model with the ansaetze 

£>/,+-,++ = A [L + L s ]- 2 (38) 

and 

£>/,+-,++ = A [L + L s + d(L + Ls)- 1 ]- 2 . (39) 

Fitting all data with the ansatz (EH} we get A = 3.20901(25), L s = 1.9140(11) and x 2 /d.o.f. 
= 1.12. Fitting all data with the ansatz (EH} we get D = 3.2071(5), L s = 1.898(4), d = 
0.20(6) and x 2 /d.o.f. = 0.72, instead. We redid these fits for Df t+ _ t++ evaluated at (3 = 
0.38772176 in order to estimate the error due to the uncertainty of (3 C . Finally, in order to 
check for the possible effect of residual corrections to scaling oc Lq w , we fitted our data with 
the ansaetze (13511361) . where we fixed the amplitude of the leading correction to c = 1.5/30. 
Note that in ref. [7| we found that the amplitudes of the leading correction are at least 
suppressed by the factor 1/30 in the Blume-Capel model at D = 0.655 compared with the 
spin- 1/2 Ising model. 

Taking these fits into account we arrive at 

A + _ - A ++ = 3.208(5) (40) 

which is consistent with the estimate ([57} obtained above. Furthermore these results are fully 
consistent with A + _ - A ++ = [0+_(O) - 0++(O)]/2 = [5.613(20) + 0.820(15)]/2 = 3.217(18) 
obtained in section VI C of ref. [26j. Our result is slightly larger than A + _ — A ++ = 
2.71(2) — 0.376(29) = 3.09(5) which the authors obtained by fitting their data for the 



thermodynamic Casimir force per area with ansatz (26) of ref. [22]]. In [28J the authors 
used different ansaetze. Eqs. (17,18,19) coincide at the critical point with our ansatz (jTJ). 
The authors argue that corrections oc Lq w are effectively taken into account by the oc L^ 1 
correction that is present in the ansatz. In figure 6 a of [28] we see that their strong symmetry 



breaking results, i.e. hi = —100 and h\ = 100 clearly deviate from ours 26]. To understand 
this discrepancy we have fitted our data for the Ising model with the ansatz ( 139} . Fitting 
all our data we get A = 3.1467(4), L s = 3.480(4), d = -5.83(5), and x 2 /d.o.f.= 76.35. 
Fitting only the data with L < 34 and assuming a statistical error that is 3 times larger 
than the one that we acctually achieved we get A = 3.136(2), L s = 3.39(2), d = —4.7(2) 
and x 2 /d.o.f.= 1.03. While x 2 /d.o.f.~ 1, this is completely incompatible with our final 
result ( E7J , which substantiates our statements on fitting with appoximate ansaetze above. 

Finally note that our results for L s of the Blume-Capel model at D = 0.655 are fully 
consistent with L s = 1.9(1) 0, L s = 2l ex = 1.92(4) and L s = 1.90(5) [29]. In section [V]] 
below, we shall assume L s = 1.91(5). 
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TABLE II. Number of measurements (stat) in our simulations of the Ising model at (3 = 0.2216546. 
For each measurement 16 sweeps with the Metropolis algorithm were performed. In these simula- 
tions L = 6Lq and L = IOLq for ++ and H — boundary conditions, respectively. 



Lq stat ++ stat H — 



6 


64. 


.0 


X 


10 8 


64. 


.0 


X 


10 7 


7 


57. 


.2 


X 


10 s 


64. 


.0 


X 


10 7 


8 


45. 


.3 


X 


10 s 


64. 


.0 


X 


10 7 


9 


47. 


.9 


X 


10 s 


64. 


.0 


X 


10 7 


10 


39. 


.4 


X 


10 s 


51. 


.5 


X 


10 7 


11 


31. 


.4 


X 


10 8 


46. 


.1 


X 


10 7 


12 


24. 


.0 


X 


10 8 


44. 


.8 


X 


10 7 


13 


15. 


.1 


X 


10 8 


37. 


.9 


X 


10 7 


14 


15. 


.5 


X 


10 s 


32. 


,4 


X 


10 7 


15 


15. 


.3 


X 


10 s 


27. 


.8 


X 


10 7 


16 


14. 


.2 


X 


10 s 


25. 


.6 


X 


10 7 


17 


10. 


.4 


X 


10 s 


21. 


.5 


X 


10 7 


18 


10. 


.9 


X 


10 8 


19. 


.5 


X 


10 7 


19 


11. 


.8 


X 


10 8 


18. 


.9 


X 


10 7 


20 


10. 


.4 


X 


10 8 


18. 


.8 


X 


10 7 


22 


10. 


.4 


X 


10 8 


28. 


.7 


X 


10 7 


24 


67. 


.1 


X 


10 7 


20. 


,4 


X 


10 7 


26 


64. 


.3 


X 


10 7 


22. 


.7 


X 


10 7 


28 


62. 


.0 


X 


10 7 


25. 


.8 


X 


10 7 


32 


62. 


.9 


X 


10 7 


22. 


.5 


X 


10 7 


36 


31. 


.5 


X 


10 7 


22. 


.7 


X 


10 7 


48 


14. 


.4 


X 


10 7 


2. 


.9 


X 


10 7 


64 


9 


.9 


X 


10 7 


2. 


.7 


X 


10 7 



B. Simulations at the critical point 

In order to compute the energy per area and the magnetisation profile at the critical 
point of the Ising model, we performed high statistics simulations at — 0.2216546, which 
was our estimate of (3 C when we started the simulations. In order to obtain the observables 
at /3 = 0.22165462, we computed the derivate of the observables with respect to /3 from 
finite differences. In table [Til we summarize the lattice sizes and the statistics of our first set 
of simulations. In a second set of simulations with H — boundary conditions we measured 
the second moment correlation length in addition. We simulated lattices of the thicknesses 
L = 24, 32, 48, 64 and 96. The number of measurements is 51.2 x 10 7 , 51.9 x 10 7 , 49.1 x 10 7 , 
34.6 x 10 7 , and 7.8 x 10 7 , respectively. Also here we performed 16 sweeps with the Metropolis 
algorithm for each measurement. For this second set of simulations L = 4Lq. For Lq = 6 
we simulated L = 12, 14, 16, 18, 20, 24, 36, 48 and 60 performing 6.4 x 10 8 measurements 
throughout. From the analysis of these runs we conclude that for H — boundary conditions, 
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at the critical point L = ALq is fully sufficient to keep deviations from the L — > oo limit at 
a negligible level. In our simulations we wrote averages over 64000 measurements on disc to 
keep the amount of data tractable. In order to estimate autocorrelation times we did a few 
additional simulation, where every measurement was stored. For example we performed 10 5 
measurements for H — boundary conditions, L = 96 and L = 384. From this run we got the 
integrated autocorrelation times r int = 3.3(2), 15.2(1.0) and 28.(3.) in units of measurements 
for the energy per area, the magnetic susceptibility and the magnetisation in the middle of 
the film. The autocorrelation times of a local algorithm grow like r oc Lq at the critical 
point, where z ~ 2. Therefore, despite the efficient multispin coding implementation of the 
Metropolis algorithm, the cluster algorithm should become more efficient starting from a 
certain thickness Lq. Since Tj ni enters into the statistical error this thickness depends to 
some extend on the observable one is interested in. For lack of human time, we did not 
systematically investigate these questions. 

C. The energy per area 

Taking eq. ( 130]) at x = and ignoring corrections to scaling we arrive at 

E ex (L ,P c ) = B + aL 2+1/u (41) 

where B = 2f r (/3 c ) and a = £~ 1/v h'(0). 

In order to compute the excess energy, we used the estimate of E bu i k ((3 c ) obtained in the 
appendix. Replacing L by £o,e// iu eq. fl4"Tl) we arrive at the ansaetze 

E ex (L , (3 C ) = B + a[L + L s + c(L + L^}' 2 ^ (42) 

and 

E ex (L , (3 C ) = B + a[L + L s + c(L + L s ) x ~ u + d(L + L s ) 1 ~ e ]~ 2+1 / V (43) 

where we set either e = 1.664 or e = 2. In our fits, B, a, c, L s and d are free parameters. 
We fixed v = 0.63002 and u = 0.832. 

First we analysed our data for H — boundary conditions. Fitting with the ansatz (T4"2l we 
get an acceptable x 2 /d.o.f. starting from L 0t7nin = 18. For L 0t7nin = 20 we get B = 7.8010(7), 
a = -15.455(7), c = 1.472(35), L s = 1.413(44) and x 2 /d.o.f.= 0.62. Using the ansatz (PJ 
we get an acceptable x 2 /d-o.f. already for Lo,mm = 7 both for e = 1.664 and e = 2. 
For example for L^ min = 8 and e = 1.664 we get B = 7.80405(23), a = -15.4946(19), 
c = 2.028(8), L s = 6.476(10), d = -0.27(4) and x 2 /d.o.f.= 0.93. Instead for L^ min = 8 and 
e = 2 we get B = 7.80279(25), a = -15.4790(20), c = 1.794(9), L s = 0.903(11),' d = 0.13(4) 
and x 2 /d.o.f.= 1.13. We see that the results depend strongly on the ansatz that is used. 
This holds in particular for the estimates of c and L s . We redid the fits using shifted values 
of the input parameters to estimate the error of our results due to the uncertainty of these 
parameters. Taking into account the results of all these fits we arrive at B = 7.803(5) and 

a 7i+ _ = -15.48(5) - 130 x {v - 0.63002) (44) 

where for a + _ we give the dependence on the value of v explicitly. The error induced by the 
uncertainty of the other input parameters is included into the number given in (). 

For ++ boundary conditions fitting with the ansatz fj4*2|) gives acceptable values of 
X 2 /d.o.f. already for Lo imin = 7. For example for Lo,mm = 8 we get B = 7.80168(22), 
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a = -10.2105(16), c = 1.462(7), L s = 1.16(7) and x 2 /d.o.f.= 1.23. Instead fitting with the 
ansatz (g3]) we get for e = 1.664 and L %min = 8 the result B = 7.8054(7), a = -10.248(6), 
c = 2.02(4), L s = 0.27(5), d = 2.3(4) and x 2 /d.o.f.= 1.02. Fixing e = 2 we get results 
that lie between those of the two fits discussed before. Also in the case of ++ boundary 
conditions we redid the fits with shifted values of the input parameters. As our final result 
we quote B = 7.804(5) and 

a I;++ = -10.23(5) - 70 x [y - 0.63002) . (45) 

Note that the results obtained for B with H — and ++ boundary conditions agree as theo- 
retically expected. 

Assuming v = 0.63002 we get h' + _(0) = -15.48(5) x 0.1962(1) 1 / - 63002 = -1.167(5) and 
h' ++ (0) = -10.23(5) x 0.1962(l) 1 /°- 63002 = -0.771(5) from the analysis of the Ising model. 
In ref. [26[ we found for the Blume-Capel model at D = 0.655 the results obc,++ = —8.04(1) 
and a,Bc,+- = —12.18(3) for ++ and H — boundary conditions, respectively. Hence h' + _(0) = 
-12.18(3) x 0.2283(l) 1 /o-630o 2 = 1.168(4) and h' ++ (0) = -8.04(1) x 0.2283(l) 1 /°- 63002 = 
—0.771(2). We see that the results obtained for the universal quantities h' + _(0) and h' ++ (0) 
are in perfect agreement. Using eq. ( !29l) we arrive at 

9' + _(0) = -0.482(2) , 0^(0) = -0.318(2) (46) 

taking into account the results obtained from both models. 

Finally we analysed the difference De,+-,++ at the critical point. The advantage of this 
quantity is that the bulk energy and the surface contributions exactly cancel. We fitted our 
data with the ansaetze 

E ex (L , (3 C ) = D a [L + L S + c(L + L,) 1 -]-^ (47) 

and 

E ex (L , C ) = D a [L + L S + c(L + + d(L + L s ) 1 ~ € \~ 2+X ^ v . (48) 

Fitting with the ansatz (|47p we get an acceptable x 2 /d.o.f. only for rather large values of 
L 0)min . For example for L , mi „ = 26 we get D a = -5.2548(19), c = 1.80(10), L s = 1.49(14) 
and x 2 /d.o.f.= 1.26. Fitting with the ansatz (I48p and e = 1.644 we get for I/o, m j n = 12 the 
results D a = -5.2570(8), c = 2.15(4), L s = 0.81(5), d = -3.56(18) and x 2 /d-c.f.= 1.04. 
Instead for e = 2 we get D a = -5.2548(7), c = 1.92(3), L s = 1.25(4) and d = -3.71(16) and 
X 2 /d.o.f.= 1.09. Also here, we redid the fits with shifted values of the input parameters. We 
arrive at the final result 

0/,+- - 0/,++ = -5.256(4) - 75 x [y - 0.63002) . (49) 
D. The magnetisation profile 

For simplicity, we shall not study the complete magnetisation profile, but we shall restrict 
ourselfs on the magnetisation in the middle of the film and the slope of the magnetisation 
in the middle of the film in case of ++ and H — boundary conditions, respectively. 

Let us first discuss the case of ++ boundary conditions. The magnetisation in the middle 
of the film at the critical point behaves as 

m-mid — C m L Q ^ . (50) 
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The amplitude C m is not universal, but one can construct universal amplitude ratios that 
combine C m with the amplitude of the bulk correlation length and the bulk magnetisation 
or the magnetic susceptibility. Here we only intend to compare our result for C m j for the 



Ising model with C mj Bc obtained previously for the Blume-Capel model at D = 0.655 [26 
To this end it is sufficient to determine the relative normalization of the magnetisation 
between these two models. To this end we compare the magnetic susceptibility of systems 



with the extension L = L\ = L 2 and 
that we computed in relation with ref. 



periodic boundary conditions in all three directions 
7[. In particular we fitted the data for the magnetic 



susceptibility at Z a /Z p = 0.5425 with the ansatz 

X = C X L 2 ~ V x (1 + cL- w ) + b (51) 
where we fixed rj = 0.03627(10) and u = 0.832(6). We arrive at 



C 

X)I - 1.2811(2) (52) 



C x,bc 

where statistical and systematical errors as well as the uncertainty of rj and u are taken into 
account. 

In order to define the magnetisation in the middle of the film for even values of the 
thickness L Q we quadratically extrapolated the magnetisations of the slice that is next to 
the middle and the one that is next to next. We fitted these data with the ansatze 

mmid = C m (L + L s + c(L + L^)-^ (53) 

or 

m mid = C m (L + L S + c(L + L s ) l ~" + d(L + L,) 1 ^)"^ . (54) 

where we fixed (3/v = (1 + rj) = 0.5018135, u = 0.832 and e = 1.664 or e = 2. In the 
following we only take into account data for even values of L . Using ansatz (I5"3"j) we get for 
Lo.min = 24 the results C m = 1.71799(18), c = 1.63(2), L s = 1.31(3) and x7d.o.f.= 0.21. 
Using ansatz ( IS"!]) we get with e = 2 an acceptable x 2 /d.o.f. already for L 0rnin = 6. For 
Lo, mi n = 8 we get the results C m = 1.71880(7), c = 1.844(8), L s = 0.922(11), 'd = 1.289(14) 
and x 2 /d-o.f.= 0.78. For e = 1.664 and L ,mm = 10 we get C m = 1.71929(11), c = 2.016(17), 
L s = 0.563(29), d = 1.172(27) and x 2 /d.o'.f.= 0.59. 

We redid these fits using shifted values of (3 C , rj and u. As final results we quote 

C mJ = 1.7187(10) + 4.8 x (rj - 0.03627) (55) 

where we give explicitly the dependence of our result on the value of rj. 



In ref. [26] we analysed m m ^ for the Blume-Capel model at D = 0.655 for thicknesses 
up to L = 32. Later [29j we added data for L = 48, 64 and 96. Taking into account also 
these data we arrive at 



We get 



C m:BC = 1.3421(8) + 2.8 x (rj - 0.03627) (56) 
C 

— — = 1.2806(16) (57) 

C m ,BC 
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which is fully consistent with eq. ( 152]) . 

In the case of H — boundary conditions, we consider the slope of the magnetisation profile 
in the middle of the film. It scales as 

S mid = C S L~ 1 -^ . (58) 

We fitted our data for the Ising model with the ansaetze 

S mid = C S (L + L S + c(L + L,) 1 -")- 1 -^ (59) 

and 

S mid = C S (L + L S + c(L + L s y-"d{L + Lsf-'Y 1 '^ (60) 

where we fixed rj = 0.03627 and u = 0.832 and e = 1.664 or e = 2. Also here we fitted 
only the data for even values of L . Fitting with the ansatz (159"]) we find small values of 
X 2 /d.o! already for L 0>rnin = 8. For L , m in = 10 we get C sJ = 7.2013(4), c = 1.4603(25), 
L s = 0.7023(31) and x 2 /d.o.f.= 0.39. Fitting with the ansatz (160]) we find that the parameter 
c vanishes within the error bars. Taking into account the error due to the uncertainty of the 
input parameters u and 77 we arrive at the 

C sJ = 7.201(3) + 19 x (77 - 0.03627) . (61) 

Fitting data obtained in relation with ref. [26[ for the Blume-Capel model we get 

C S>BC = 5.625(10) + 10 x (77 - 0.03627) . (62) 

We get 

C 

= 1.280(3) (63) 

^s,BC 

which is fully consistent with eq. ( 152"]) . 



E. The correlation length 

Finally we discuss the second moment correlation length of films with +— boundary 
conditions at the critical point. Our numerical results are summarized in table IIHI Since 
here we generated less data than for the quantities discussed above we abstain from fitting 



the data for the correlation length. In ref. |26j we found ^2nd = 0.2115(8)(L + L s ). Based 
on this result we define 

L ,ef f = £W0.2115(8) . (64) 

In the third column of table [Till we quote £o,e// — L . In [] we give the error due to the 
uncertainty of the amplitude of the correlation length of the film. For comparison we give 
analogous results derived from the difference of free energies _D/ >+ _ ;++ , the difference of 
energies D Ej+ _ j++ , the magnetisation in the middle of the film for ++ boundary conditions 
and the slope of the magnetisation in the middle of the film for H — boundary conditions. 

We see that the values of -^o,e// — Lq computed from different observables are of a similar 
size. However the differences are considerably larger than the sum of the errors. Therefore 
it is quite clear that -^ , e // — L is not exactly the same for all quantities. 
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TABLE III. In the second column we give the second moment correlation length obtained from 
simulations of lattices with L = 4Lo for H — boundary conditions at the critical point of the Ising 
model. In the third column we give L ex = Lq^jj — Lq. For the definition of Lo,eff see the text. 

In the fourth, fifth, sixth, and seventh column we give L ex = -^o,e// — Lq derived from Df _| 4.+, 

De,+-,++, m mi d, and S mid , respectively. 



Lo (2nd L ex , (2nd 


L e x, Df^ — ^ 


+ L ex , De,-\ — ,++ 


Lex, n^mid 


Lex, ^mid 


24 5.6881(24) 2.89[10] 


3.51[2] 


4.61[10] 


4.14[3] 


3.20[1] 


32 7.4025(42) 3.00[13] 


3.64[3] 


4.76[12] 


4.27[4] 


3.32[1] 


48 10.807(10) 3.10[19] 


3.83[4] 


4.99[16] 


4.49 [6] 


3.51[1] 


64 14.204(20) 3.16[25] 


3.97[5] 


5.16[20] 


4.65[8] 


3. 64 [2] 


96 20.99(10) 3.2[4] 








3.85[3] 



V. THERMODYNAMIC CASIMIR FORCE AND THE TRANSFERMATRIX 



First let us briefly recall the discussion given in section IV of ref. [26(. The partition func- 
tion of a system with fixed boundary conditions can be expressed in terms of the eigenvalues 
of the transfermatrix and the overlap of the eigenvectors with the boundary states. Let us 
consider a lattice of the size Lq x L 2 , where L is large compared with the bulk correlation 
length but still finite. We consider the transfermatrix T that acts on vectors that are build 
on the configurations living on L 2 slices. We denote the eigenvalues of T by X a and the 
corresponding eigenvector by |a), where a = 0, 1, 2, a max . The eigenvalues are ordered 
such that X a > \p for a < (3. In particular A is the largest eigenvalue. The partition 
function of the system with fixed boundaries is given by 



'h,b 2 



J2^ a (bi\a)(b 2 \a) , (65) 



where / = Lq + 1 for our definition of the thickness Lq. The boundary states are either 
+ or — here. It follows that 

L 2 d 

;Fcasimir = 7^ [ m ^61,62 — ^ m ^0] 



£ a ln(A Q /A ) (XJXq) 1 {h 


a) (62 


a) 




a) (b 2 


a) 



Y, a m a exp(-mj) {h 


a) (6 2 


a) 




a 


)(b2 


a 


) 



KbT ol 

(66) 

where l/£ a — m a = — ln(A a /Ao) are inverse correlation lengths. In the high temperature 
phase for ^ = ^ C L the force is dominated by the contribution from a = 1. Hence 

9W * ^F C0 . jmir « - m Vexp(- mi ) j^fgfg ■ (67) 
The finite size scaling behaviour of the thermodynamic Casimir force implies that 

mL(b\0) v ; 
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has a finite scaling limit. The state |0) is symmetric under the global transformation s x — > 
—s x for all x in a slice, while |1) is anti-symmetric and therefore C = C + = —C-. Hence 

e ++ (ml) = -0 + _(mZ) = -C 2 m 3 Z 3 exp(-m/) (69) 

for sufficiently large values of ml. Since x = t[l/£o] 1/u ^ {ml) 1 '" it follows 

e ++ (x) = -6+-(x) ~ -C 2 x 3u exp{-x u ) (70) 

for sufficiently large values of x. 



A. C and the magnetisation profile 

In the following we shall discuss how the overlap amplitude C 2 can be computed from 
the magnetisation profile of a semi-infinite system with + boundary conditions and the cor- 
relation function of slice magnetisations. In terms of the transfermatrix, the magnetisation 
at position xq in a film of thickness Lq is given by 

(M x = ( > s X0yXuX2 ) = ^ 71 



\X1,X2 



In the basis of slice configurations, M is a diagonal matrix, where the elements give the 
magnetisation of the corresponding configuration. For I ^> £ and £ 2 <C x <C / eq. (JTTJ) 
reduces to 

(M(xo)) = ArA^°(& 1 |l)(l|M|0)(0|& 2 ) = ^ ( Ai N - 



A(,(&i|0><0|& 2 > (6i|0) xI|M|0) Uo 
= mLC 6l (1|M|0) exp(-ms ) (72) 

The quantity Om = (1\M\0)/L is finite in the limit L — > oo, since (M(xq)) / L 2 is finite in 
this limit. 

The slice-slice correlation function for a lattice of linear size Lq and periodic boundary 
conditions is given by 



1 ,w, ™ , u 1 E^</^i«>HM|/3>A^°- 



■r 



G(r) = -(M(xo)M(x + r)> = (73) 

(74) 

Since M is antisymmetric under s x — >■ — s x for all x in the slice, (0|M|0) vanishes. For 
£2 < x < L we get 

G(r) = — (0|M|l)(l|M|0)exp(-mr) = O^exp(-mr) (75) 

Taking into account the periodicity of the lattice we arrive at 

_ n2 exp(-mr) + exp(-m(L - r)) 
G(rj " ° M l + exp(-mL ) Ubj 

which we shall use in our numerical analysis below. 
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B. Numerical implementation 



In order to compute G(r) we simulated lattices with L = L\ = L 2 = L and periodic 
boundary conditions. In the case of the Blume-Capel model we simulated the mode l by 



using a hybrid 35] of the local heat-bath algorithm and the single-cluster algorithm |34|. 
In the case of the Ising model we only used the single-cluster algorithm. We measured 
the correlation function G(r) by using its cluster-improved estimator. In order to keep 
deviations from the thermodynamic limit negligible we we chose L > 10£ throughout. For 
a discussion of this point see section III or ref. [HI]. In order to compute £ and 2 M from 
eq. (176|) we took the correlation function at the distance r and r + 1. For eq. ( I75|) one gets 
£ = 1/ hi(G(r)/G(r — 1)) and 2 M = G(r) exp(r/£). For eq. ( I76|) we solved the system of two 
equations numerically. We computed the statistical errors of £ and 2 M and their covariance 
by using the Jackknife method. We checked which distance r is needed to keep corrections 
due to eigenstates of the transfer matrix with a > 1 negligible. As a result, we took r 2£ 
throughout. 

In the case of the Blume-Capel model at D = 0.655 we simulated at 11 values of (3 
between (3 = 0.34 where £ = 1.50420(13) and (3 = 0.3872 where £ = 26.7102(16). For 
(3 = 0.3872 we performed about 10 7 update cycles. Each cycle consists of two sweeps of the 
local heat-bath algorithm and 10 4 single-cluster updates. Note that the average cluster size 
at (3 = 0.3872 is 1645.58(17), and hence the lattice of the size 270 3 is covered on average 
0.84 times by these 10 4 clusters. The simulation at = 0.3872 took about 13 month of 
CPU-time on a single core of a Quad-Core AMD Opteron(tm) Processor 2378 running at 
2.4 GHz. In the case of the Ising model, we simulated at 59 values of ft between = 0.125 
where £ = 0.667308(53) and (3 = 0.2208 where £ = 16.6711(12). 

Next we analysed the magnetisation profile of films with ++ boundary conditions. Also 
here we required that Li > 10£. When possible, we used the results obtained from the 
simulations that we performed to compute the thermodynamic Casimir force. For values of 
(3 where this is not the case, we performed extra simulations using the cluster algorithm. 
Taking 2 M and £ obtained above from the simulations of the lattices with periodic boundary 
conditions as input one gets an estimate of C(£) from eq. (!72|) for each distance Xq from the 
boundary. Throughout we took our final result from Xq ~ 3£. 

In figure [T] we plot our results for C(£) as a function of m = l/£ for the Ising model 
and the Blume-Capel model at D = 0.655. Note that the error bars are much smaller than 
the size of the symbols. For example for the Blume-Capel model at (3 = 0.3872 we obtain 
C(f) = 1.2241(4) and for the Ising model at (3 = 0.2208 we get C(f) = 1.1500(3). The data 
for the Blume-Capel model essentially fall on a straight line, confirming that corrections 
oc are eliminated and those oc caused by the boundaries dominate. In contrast, 
for the Ising model we see a clear bending of the curve. It is conceivable that in the limit 
£ — > oo the two curves converge to a unique value. 

In order to substantiate these qualitative observations we we fitted our data with the 
ansaetze 

C(0 = Cexp(-c/0 (77) 

and 

C(0 = Cexp(-c/0 + ar e (78) 

where C, c and a are the parameters of the fit. First we analysed our data for the Blume- 
Capel model. Fitting with the ansatz (1771) we get x 2 /d.o.f.= 0.67, for fitting all data except 
the smallest value of 0. The results for the parameters of the fit are C = 1.24568(21) 
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FIG. 1. The amplitude C(£) for the Ising model and the Blume-Capel model at D = 0.655 as a 
function of l/£. 



and c = 0.4572(7). Next we fitted all data with the ansatz ( J78|) . Fixing e = 0.832, we 
get C = 1.2462(5), c = 0.442(9), a = -0.017(11) and x 2 /d.o.f.= 1.06. For e = 2 we get 
C = 1.24588(27), c = 0.4591(15), a = 0.0043(23) and x 2 /d.o.f.= 0.64. As our final estimate 
we give 

(79) 



C = 1.2459(7) 

where the error-bar covers the results of the three fits given above. The estimate C 
given in 
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1.5(1) 
1.552(2) 
29l Note 



is consistent with, but much less precise than our present estimate C 2 
Note that the result c « 0.46 is fully consistent with l ex = 0.96(2) obtained in 
that for our definition of the thickness one expects c = l ex — 1/2. 

Next we fitted our data for the Ising model with the ansatz flTSJ) using e = 0.832. Fitting 
all data with (3 > 0.202 we get C = 1.24653(23), a = -1.3750(29), c = -0.479(2) and 
X 2 /d.o.f.= 1.17. Taking into account smaller values of /3, x 2 /d.o.f. rapidly increases. We 
redid the fit using e = 0.826 and we also fitted with ansaetze that include subleading 
corrections. Taking into account the results of these fits we arrive at C — 1.247(3), which 
is fully consistent with the result ( I79p that we obtained from the data for the Blume-Capel 
model. 

We performed a similar study to determine the behaviour of the thermodynamic Casimir 
force for ++ boundary conditions for x — > —00 in the low temperature phase. However here 
we can not reach the same precision as above, since there is no efficient improved estimator 
for the correlation function in the low temperature phase, and secondly contributions due 
to subleading states of the transfermatix are more important than in the high temperature 
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phase. In the case of the Blume-Capel model we have computed C for 16 values of (3 in the 
range from (3 = 0.39 where £ = 5.584(40) up to f3 = 0.405 where £ = 1.5697(49). In the case 
of the Ising model in the range from = 0.223 where £ = 6.6028(20) up to (3 = 0.227 where 
£ = 2.7321(42). 

Analysing the data for the Blume-Capel model, fixing c = 0.46(2) we arrive at C — 
0.428(10) and hence C 2 = 0.183(9) which is consistent with but more precise than C 2 = 



0.20(5) given in [26] . Analysing the data for the Ising model, we get a consistent result. 



C. The correction function 



Plugging in C 2 {t) = C 2 {1 + a c t e ) and £ = 6r"(l 
H — boundary conditions 



dfex 

0L 



Lv 3 C 2 ^t u exp I -^t 



U 



L 



£o 



1 + {a c + 



d£t 6 ) into eq. (169]) we get, e.g. for 
a \ f + (t 29 ) 



with 



L 3 6(x) x [1 + bq(x)L " + O(L 

bq( x ) = $(a c + [x u - 3}a^)x e . 



which is not consistent with 



bq(x) 



-ex 



(80) 

(81) 
(82) 



that one derives by plugging eq. ( 170|) into eq. ( flOl . In figure [2] we plot q(x)x v as a function 
of x. To this end, we take the numerical values £o — 0.1962, = —0.32, eq. (lAllj) . and 
a c = 2 x 0.1962- - 832 x (-1.375/1.247) = -8.55. It turns out that the curve is very flat in 
the range of x we are interested in. Also the value is rather close to the values of c that we 
obtained from the analysis of data directly at the critical point. 



VI. THE SCALING FUNCTION OF THE THERMODYNAMIC CASIMIR FORCE 
FOR ++ AND +- BOUNDARY CONDITIONS 

We computed the thermodynamic Casimir force using the method discussed by Hucht 



2l[. Starting from the energy per area we computed 

AE ex (L , (3) = [E(L + d/2, (3) - E(L - d/2, (3)}/d - E bulk {(3) . (83) 
The value of the energy density of the bulk system E bu i k ((3) is obtained from an anal ysis 



of the high temperature series given in [39] and the low temperature series given in [40 



combined with Monte Carlo simulations. For details see Appendix IA31 

In order to obtain Af ex we numerically integrated AE ex using the trapezoidal rule: 



n— 1 

- Af ex ((3 n ) w -Af ex {(3 ) + 2^<+i ~ l&EeM+i) + A^(A)] (84) 



i=0 



where are the values of (3 we simulated at. They are ordered such that (3 i+ i > for all i. In 
previous work (3 had been chosen such that AE ex ((3 ) w and therefore also Af ex ((3 ) ~ 0. 
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FIG. 2. We plot q(x)x " as a function of the scaling variable x in the range that is relevant for 
our problem. For the definition of q(x) and a discussion see the text. 



Here, instead we chose a somewhat larger value of /3q such that the approximation discussed 
in the previous section is still valid. In particular, we set 

Afm ^ C 2 (A)) exp[-(L + 1 + d/2)/g(fl))] - exp[-(L + 1 - d/2)/g(A,)] 
fM) ~ M d (85) 

where we have the + sign for ++ boundary conditions and the — sign for H — boundary 
conditions. By comparing results obtained with different choices of flo we found that the 
approximation f l8"5"j) is accurate at the level of our statistical error up to L /£(A)) ~ 8- To 
be on the safe side, we have used Z/ /£(/3o) > 10 in the following. 

We simulated the Ising model with ++ boundary conditions for the thicknesses Lq = 8, 
9, 14, 15, 16, 17, 18, 19, 32, 34, 64, and 68. Using the resulting data we computed the 
thermodynamic Casimir force for the thicknesses Lq = 8.5 and Lq = 16.5 using the difference 
d = 1. In order to check for the effect of using a finite difference to compute d/dLo we have 
redone the calculation for Lq = 16.5 using d = 3 and 5 in addition to 1. We conclude that 
d/Lo rs 0.06 is sufficient at the level of our accuracy. Therefore for Lq = 33 and Lq = 66 we 
used d = 2 and d = 4, respectively. Throughout we used L > 5L , which is clearly sufficient 



to neglect deviations from the limit L — > oo; See ref. [26|. We chose (3q = 0.15, 0.19, 0.21 
and 0.218 for L = 8.5, 16.5, 33 and 66, respectively. We simulated at 163, 122, 117 and 41 
values of (3 for these thicknesses, respectively. Note that in the case of L = 66 we simulated 
only up to 0c, since these simulation are rather expensive. 

For Lq = 16 and 17 we performed 6.4 x 10 8 measurements for each value of ft that we 
simulated at. For each measurement we performed 16 sweeps with the Metropolis algorithm. 
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FIG. 3. We plot 0^ , — 9 ++ and the approximation (|85|) . The data are taken for the Blume-Capel 

model at D = 0.655 and the finite difference is computed from Lq = 32 and Lq = 34. 



In total these simulations took about 8 years of CPU time on one core of a Quad-Core 
AMD Opteron(tm) Processor 2378 running at 2.4 GHz. For Lq = 15 and 18 we performed 
1.3 x 10 8 measurements and for Lq = 14 and 19 only 6.4 x 10 7 measurements. For L = 32 
we performed between 2.6 x 10 7 and 6.4 x 10 7 measurements and for Lq = 34 we measured 
2.6 x 10 7 or 3.2 x 10 7 times for each value of /?. These simulations took about 5 years of 
CPU time on one core of a Quad-Core AMD Opteron(tm) Processor 2378 running at 2.4 
GHz. For Lq = 64 and 68 we performed 6.4 x 10 6 measurements for each value of /3. In 
total these simulations took about 2.5 years of CPU time on one core of a Quad-Core AMD 
Opteron(tm) Processor 2378 running at 2.4 GHz. 



We improved the numerical results obtained in ref. [26] for the Blume-Capel model. To 
this end, we simulated at additional values of 0. This way both the statistical error of our 
result as well as the systematical error of the numerical integration are reduced. In ref. 26| 



we did simulate the thicknesses L = 8, 9, 16, 17, 32 and 33. Here we simulated L = 34 in 
addition. 

In figure [3] we plot 6> + _, — 9 ++ and the approximation (185]) computed by using the data 
obtained for the Blume-Capel model at D = 0.655 for L = 33 and d = 2. As discussed 
at the end of section HV A\ we used the value L s = 1.91 to compute the effective thickness 
Lo,eff — Lq + L s . The deviation of 6* + _ and — 9 ++ from the approximation (185]) is smaller 
than 5% for x ^ 16 and smaller than 1% for x ^ 22.5. The average (0 + - — 9 ++ )/2 deviates 
from the approximation fl85l) by less than 5% for x 8.6 and by less than 1% for x ^ 12.7. 

Next we extracted the value and the location of the minimum of A/+ + . In the case of 
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-20 1/v 20 40 

FIG. 4. We plot A/Lq as a function of ^Lo.e/z/Co] 1 ^ for ++ boundary conditions. The thick 
lines give the result obtained for the Blume-Capel model at D = 0.655 and the two thicknesses 
Lq = 16.5 and Lq = 33. In the case of the Blume-Capel model we have used Lo >e ff = Lq + 1.91 
as effective thickness of the film. Our results for the Ising model are given by thin lines. In the 
case of the Ising model we have used the effective thicknesses £o,e// = 19.712, £o,e// = 36.509 and 
Lo,eff = 69.936, for Lq = 16.5, Lq = 33 and Lq = 66, respectively. These effective thicknesses are 
chosen such that at the minima the curves fall on top of the one for the Blume-Capel model and 
Lq = 33. At the resolution of the plot, all 5 curves fall on top of each other almost everywhere. 
Only for 20 ;$ x J$ 40 the curve for the Ising model and Lq = 16.5 can be distinguished from the 
other four. 



the Blume-Capel model we get (3 min = 0.382185(15) and Af ++ . min = -0.0002808(6) for 
Lq = 16.5 and (3 min = 0.385716(6) and Af ++tTnin = -0.00004117(5) for Lq = 33. This 
corresponds to t min [L 0je ff/£o] 1/l/ = 5.88(5) and Lg e// A/ ++)min = -1.752(18) for L = 16.5 
and t min [LQ :eff /^ = 5.88(4) and L^ e// A/ ++iTOi n = -1.752(10) for L = 33. The quoted 
error-bars include the error of (3 m i n , Af ++tmin and errors induced by the uncertainties of L s , 
£o, v and (3 C . The values obtained from Lq = 16.5 and Lq = 33 agree nicely. Our results are 
also consistent with those of ref. [26j]: x min = 5.82(10) and 9 ++ ^ min = —1.76(3). Our results 
obtained for the Ising model are summarized in table IIVI Here we have computed £o,e// 
by requiring Ll e jj-Af min = —1.75169... which is our estimate obtained for the Blume-Capel 
model and L = 33. We see that the values of L Q e ff are similar to those obtained in section 
|IV]from the analysis of the free energy differences at the critical point. In the last column 
we give t m m[-^o,e// Hq} 1 ^ u using these values of L 0>e ff. We see that these estimates of x m i n 
are essentially consistent with that obtained above from the analysis of the Blume-Capel 
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TABLE IV. Results for the minimum of 9 ++ obtained for the Ising model. 



L - d/2 L 


+ d/2 






L 0,eff 


tmin[Lo,eff /{.o] 1 ^ 


8 


9 


0.2123025(16) 


-1.1605(1) xlO" 3 


11.471 


5.96(1) 


14 


19 


0.2176215(5) 


-2.347(1) xlO~ 4 






15 


18 


0.2176744(19) 


-2.306(1) xlO" 4 






16 


17 


0.2176975(30) 


-2.2869(15) xlO~ 4 


19.712 


5.96(1) 


32 


34 


0.2201704(30) 


-3.5996(26) xl0~ 5 


36.509 


5.94(2) 


64 


68 


0.2211284(25) 


-5.121(18) xlO" 6 


69.936 


5.91(3) 



model. 

For L = 16.5 we have checked the effect of the discretization error on the position and 
the value of the minimum. The error behaves as e = ad 2 + 0(c? 4 ). The results obtained for 
d = 1, 3 and 5 are consistent with a quadratic behaviour. For d = 1, the relative error is 
about one permille for both Af min and t m i n . 

In figure H] we plot our numerical results for the scaling function 9 ++ which are given 
by Ll e jjAf ++ as a function of tfLce/z/^o] 1 ^ where v = 0.63002 is set. In the case of the 
Blume-Capel model we use £o,e// = Lq + 1.91 as effective thickness of the film. We give our 
results for Lq = 16.5 and 33. For the Ising model we take the effective thicknesses given in 
the sixth column of table [TV] We plot our results for Lq = 16.5, d = 1, Lq — 33 and Lq = 66. 
The error bars are too small to be visible in the plot. At the resolution of the plot, all 5 
curves fall on top of each other almost everywhere. Only for 20 J$ x J$ 40 the curve for the 
Ising model and Lq = 16.5 can be distinguished from the other four. 

Next we discuss our numerical results for the scaling function 6 + _. In figure [5] we plot 
A/Lq ^ as a function of £[L ,e///6o] 1 ^' f° r the Blume-Capel model at the thicknesses Lq = 
16.5 and 33 and the Ising model at Lq = 16.5, 33 and 66. In the case of the Blume-Capel 
model we use £o,e// = L + L s with L s = 1.91. For the Ising model we take the same values 
for L 0iE ff as above for ++ boundary conditions. 

We find that the different curves fall quite nicely on top of each other. In the neighbour- 
hood of the maximum the curve for the Ising model at L = 16.5 lies slightly below the 
other ones and for x J$ —30 the curves slightly fork. The discrepancies discussed for ++ 
boundary conditions in the range 20 J$ x J$ 40 are also present for H — boundary conditions. 
They can not be seen in figure [5] since the range of values for H — boundary conditions is 
larger than that for ++ boundary conditions. 

In table |V] we have summarized results for the maximum of 6+-. Using L s = 1.91 in 
the case of the Blume-Capel model we get nicely consistent results for x max and mr , x 
from the two thicknesses L = 16.5 and L = 33. These results improve those of ref. [261 ]: 
x+-,max = —5.17(7) and 9 + _ tmax = 6.56(10). In the case of the Ising model we use the values 
of £o,e// obtained above for films with ++ boundary conditions. The resulting estimates for 

x m ax and 9^ >max are close to those obtained from the Blume-Capel model. In particular 

the results obtained for Lq = 33 are closer to the Blume-Capel ones than those obtained for 
L Q = 16.5. 

We conclude that our numerical results for the scaling functions of the thermodynamic 
Casimir force for ++ and H — boundary conditions are fully consistent with the universality 
hypothesis. Furthermore our ansatz ()8]) provides a good approximation of the universal 
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FIG. 5. We plot A/Lpgyj as a function of t[-ko,e///£o] 1 / 1 ' f° r H — boundary conditions. The thick 
lines give the result obtained for the Blume-Capel model at D = 0.655 and the two thicknesses 
Lq = 16.5 and Lq = 33. In the case of the Blume-Capel model we have used A),e// = Lq + 1.91 
as effective thickness of the film. Our results for the Ising model are given by thin lines. In the 
case of the Ising model we have used the effective thicknesses £o,e// = 19.712, £o,e// = 36.509 and 
Lo,eff = 69.936, for Lq = 16.5, Lq = 33 and Lq = 66, respectively. These values are taken from the 
analysis of ++ boundary conditions above. At the resolution of the plot, all 5 curves fall on top 
of each other almost everywhere. Near the maximum the curve for the Ising model and Lq = 16.5 
stays slightly below the other ones. For x J$ —30 the curves slightly fork. Note that in this range 
the difference between the Blume-Capel results for Lq = 16.5 and Lq = 33 is of a similar size as 
the one between the Ising results for Lq = 16.5 and Lq = 33 and between Blume-Capel and Ising. 

correction function. 

VII. SUMMARY AND CONCLUSIONS 

We studied the spin-1/2 Ising model and the improved Blume-Capel model on the simple 
cubic lattice with film geometry. In particular we considered strongly symmetry breaking 
++ and H — boundary conditions. We focused on the thermodynamic Casimir force. At 
the critical point we studied the behaviour of the free energy per area, the energy per area, 
the magnetisation profile and the second moment correlation length of the film. The main 
subject of the present work are corrections to scaling. Previously it has been demonstrated 
at the example of improved models that corrections oc L 1 that are due to the boundaries 
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TABLE V. Results for the maximum of 0_| obtained for Blume Capel (BC) model and the Ising 

(I) model. In the second and third column we give the thicknesses that have been considered. In 
the fourth column we give the value of A/ at the maximum and in the fifth column we give the 

location of the m&ximum. In th.6 sixth 3>nd SGVGnth. column wg givG Gstim&tGS ofQ-\ } ma% &nd x max 

dGrivGd from thGSG tgsuHs. 



Model L 


- d/2 L 


+ d/2 


Pmax ^fmax J-'q e f f^fmax 


t-tnax [-ko,e///£o] 


BC 


16 


17 


0.39257(3) 0.0010501(7) 6.552(5)[54] 


-5.15(3) [3] 


BC 


32 


34 


0.389474(5) 0.00015426(5) 6.563(2)[28] 


-5.139(15)[15] 


I 


16 


17 


0.224948(4) 0.00085044(30) 6.514(2) 


-4.959(6) 


I 


32 


34 


0.2229119(3) 0.000134650(35) 6.552(2) 


-5.035(12) 



can be expressed by an effective thickness £o,e// = L + L s , where L s is the same for all 
quantities. Note that L s depends on the model and in particular on the details of the 
boundary conditions. Here we have probed the hypothesis that the leading bulk corrections 
can be expressed in an analogous way: 



J 0,eff 



L + L s + c(L + (86) 



Fitting various quantities at the critical point of the Ising model we find similar, but likely 
not identical values of the amplitude c. Also the study of the thermodynamic Casimir force 
for large values of the scaling variable x shows that eq. (I8"6"j) can not be exact. Nethertheless it 
turns out to be a surprisingly good approximation in the range of x that is of experimental 
interest. In section IVII we investigate the thermodynamic Casimir force for ++ and H — 
boundary conditions. We find for AfL^^f plotted as a function of f[Lo >e ///Co] a good 
collapse of the data for both the spin- 1/2 Ising model and the Blume-Capel model. In the 
case of the Blume-Capel model we have used A),e// — L + L s with L s = 1.91(5). We have 
demonstrated that in the case of the spin- 1/2 Ising model approximately the same £o,e// can 
be used for ++ and H — boundary conditions. The values of £o,e// that we have obtained 
in section [VTl for L = 16.5, 33 and 66 are similar to those obtained from the analysis of 
Df !+ _ t++ in section HV Al The estimates of L s and c obtained from this analysis are highly 

anti-correlated. From the analysis of Df ^ j++ we get L s = 0.9 and c = 1.5 as central 

estimates. The range of possible values is given by L s = 1.1, c = 1.4 one side and L s = 0.8, 
c = 1.6 at the other. Note that the value of L= depends on the definition of the thickness. 
In particular, when comparing with refs. [22|, SO (VGMD) one should take into account 
that L y GMD = L 0ours + 2 and hence L s y GMD = L Sfiurs — 2. Since the correction function 
q(x) is universal, also for experimental data or data obtained from the numerical study of 
other models an effective thickness (|86[) should parametrize leading corrections quite well. 
Note again that L s should depend on the microscopic details of the system. In the case of 
the amplitude c universal ratios can be constructed. For example 

-8(2) (87) 



where we have used the numerical values of and £o obtained in the Appendix. In the 
introduction we have argued that eq. (jHJ) provides a good approximation for corrections to 
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scaling since fluctuations are strongly suppressed near the boundaries of the film. Therefore 
eq. OH]) should not work for periodic and anti-periodic boundary conditions. 

Furthermore we have improved the numerical accuracy of the estimates of the universal 
scaling functions 6 ++ and 6* + _: 

Writing the partition function in terms of eigenvalues and eigenstates of the transfer 
matrix and boundary states one finds for large values of x 

++ {x) = -0+_(x) = -C 2 x 3u exp{-x u ) . (88) 

Here we have demonstrated how C 2 can be accurately computed by analysing the magneti- 
sation profile of films and the bulk correlation function. We find 

C 2 = 1.552(2) . (89) 

This result can be compared with C 2 = 1.5(1) obtained in ref. 26| . 



At the critical point we find by studying the difference of free energies per area 

A + _ - A ++ = [0 + _(O) - 0++{O)]/2 = 3.204(5) (90) 

where we have averaged the results obtained from the analysis of the spin- 1/2 Ising and the 
improved Blume-Capel model. For the slope of the scaling function at the critical point we 
find 

0+_(O) = -0.482(2) , 0++(O) = -0.318(2) . (91) 

The minimum of 6 ++ is located at x m i n = 5.88(4) and takes the value # ++jm j n = 
— 1.752(10). For the maximum of we get i™ ai = —5.14(3) and + _ jTOaa: = 6.56(3). 
The reduction of the error compared with ref. 2g| is mainly due to the fact that here we 



assume L s = 1.91(5) instead of L s = 1.9(1) as in ref. [26]. 
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Appendix A: Numerical results for the spin- 1/2 Ising bulk system 
1. The critical point 

We extended the study of ref. [tJ by simulating the Ising model on the simple cubic 
lattice on a system of the size L 3 with L = 400 and periodic boundary conditions in all 
three directions at /3 — 0.2216546. As in ref. 0] we simulated the model by using a hybrid 
of the local Metropolis algorithm, the single cluster algorithm j34[ and the wall cluster 
algorithm [4l[. For details see section IV of ref. 0. We performed 2.3 x 10 7 measurements. 
In total these simulations took about 4 years of CPU time on a single core of a Quad-Core 
AMD Opteron(tm) Processor 2378 running at 2.4 GHz. In the first step of the analysis we 
determined /3 C by analysing the behaviour of the renormalization group invariant quantities 
Z a /Z p , ^2nd/L, U4 and Uq. For the definition of these quantities see section II of ref. [7]. We 
fitted our data for the Ising model with the ansatz 

R(/3 C , L) = R* + aL- w + bL~ 2 (Al) 
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where R denotes one of the renormalization group invariant quantities. Performing these 
fits, we used the results for R* given in table V of ref. 0] as input. Furthermore, we fixed 
u = 0.832. We get acceptable x 2 /d.o.f. for fits with L min > 16. The statistical error of 
(3 C increases only slowly with increasing L min . Based on fits with L min > 24 for Z a /Z p 
and ^2nd/ L we arrive at /3 C = 0.22165462(2). Instead, analysing U± and U§ we arrive at 
(3 C = 0.22165463(2). In ref. 42| the authors computed the Binder cumulant U4 on lattices 
of a linear size up to L = 1536. Fitting their data, taking the value for V\ = 1.6036(1) Q 
as input, I arrive at (3 C = 0.221654615(10). In this work we shall use 

p c = 0.22165462(2) . (A2) 

This estimate can be compared e.g. with the previous estimates /3 C = 0.22165463(8) obtained 
in ref. using a linear lattice size up to L = 96 and /3 C = 0.22165455(3) given in table X 
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At the critical point the energy density behaves as 

Ebuik(L) = E ns + aL 3 - 1 /" x (1 + bL~" + ...) (A3) 

performing various fits based on eq. ( 1A3j) . using the data of ref. ([jj]) and our new result for 
L = 400, we arrive at 

E ns = 0.9906065(15) + 85 x - 0.22165462) (A4) 

The specific heat behaves as 

C bulk {L) = C ns + aL 3 ~ 2 ^ x (1 + bL~ u + ...) (A5) 

performing various fits based on eq. (IA3j) . using the data of ref. (0]) and our new result for 
L = 400, we arrive at 

C ns = -29.1(3) - 7700000 x (p c - 0.22165462) - 3300 x (v - 0.63002) . (A6) 



2. Amplitudes and amplitude ratios 

We simulated the three-dimensional Ising model for a large number of /3-values in the 
high and the low temperature phase on L 3 lattices with periodic boundary conditions in all 
three directions. We have chosen the linear lattice size such that L > 10C,2nd{P) in order to 
keep deviations from the thermodynamic limit sufficiently small to be ignored in the analysis 
of the data. For the precise definition of the observables see section II of j3l[ . In the high 
temperature phase we simulated at 68 values of /3 in the range 0.125 < /3 < 0.2213. To 
give the reader an impression of the quality of the data, we give the results for the 5 largest 
values of (3 in table IVI1 Analogous results for the low temperature phase are given in table 

First we fitted our data for the second moment correlation length in the high temperature 
phase using the ansaetze 

6nd = 6^0,+*"" X ( X + a ^+ te ) ( A7 ) 
6nd = i2nd,0,+t- U X (1 + CL^+t 9 + bt) (A8) 

and 

&nd = 6ndA+t~" x (1 + a^ + t e + bt + ct 2v ) (A9) 
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TABLE VI. The second moment correlation length (2nd: the magnetic susceptibility x an( f the 
energy density Ebulk for the five largest values of the inverse temperature f3 that we have simulated 
in the high temperature phase of the Ising model. We have simulated 1? systems with periodic 
boundary conditions in all three directions. 

(3 L (2nd X Ebulk 

0.2206 200 14.57699(31) 831.162(32) 0.96369936(90) 
0.2207 200 15.5321(10) 940.79(11) 0.9656874(29) 
0.2208 200 16.6644(11) 1079.27(14) 0.9677195(31) 
0.2210 300 19.73548(63) 1501.960(86) 0.97198710(87) 
0.2213 400 29.1058(11) 3212.44(23) 0.97909806(69) 



TABLE VII. The second moment correlation length (2nd, the magnetic susceptibility x, the mag- 
netisation m and the energy density Ebulk for the five smallest values of the inverse temperature 
/3 that we have simulated in the low temperature phase of the Ising model. We have simulated L 3 
systems with periodic boundary conditions in all three directions. 



P L ( 2 nd X rn Ebulk 

0.2219 300 18.930(40) 1058.49(66) 0.1815607(39) 1.0126483(10) 
0.2220 200 15.294(24) 690.78(38) 0.2027298(54) 1.0200656(17) 
0.2221 200 12.976(28) 501.95(30) 0.2200006(48) 1.0271260(16) 
0.2222 170 11.418(17) 389.43(17) 0.2347800(43) 1.0339257(16) 
0.2223 170 10.278(13) 315.26(12) 0.2477779(38) 1.0405068(16) 

where t = f3 c - (3. We fixed & = 0.22165462, v = 0.63002 and u = 0.832. Based on a large 
number of fits using these ansaetze we conclude 

6nd,o,+ = 0.1962(1) +540 x (&-0.22165462) - 1 .8 x (z/-0.63002) -0.002 x (w-0.832) (A10) 

and 

a e>+ = -0.32(3)-120000x(/3 c -0.22165462) + 130x(z/-0.63002)-l.lx(w-0.832) . (All) 

Our result is in nice agreement with that of ref. [44| obtained by analysing the high tem- 
perature series of (2nd- In table VII of [44| the authors quote £o,+ = 0.5070(5) for the 
definition t = (f3 c — (3)/(3 c of the reduced temperature. Converting to our convention one 
gets £ ,+ = 0.5070(5) x 0.22165462 63002 = 0.1962(2). 

In a similar way we analysed the second moment correlation length in the low temperature 
phase and the magnetic susceptibility in both phases. Let us summarize the final results: 

6nd,o- = 0.1015(2)-200x(/3 c -0.22165462)-0.9x(z/-0.63002)-0.001x(u;-0.832) (A12) 

and 

= -0.55(15) + 70000 x [fi c -0.22165462) + 100 x(v- 0.63002)- 2.2 x(u- 0.832) (A13) 
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Using the results (1A10I) and (1A12I) we get for the universal ratio Z2nd.o.+ / £.2 nd,() .- = 1.933(5), 
which is fully consistent with £,2nd,o,+/£,2nd,o- = 1.939(5) obtained in ref. 3JJ by analysing 
Monte Carlo data obtained for the Blume-Capel model at D = 0.655. 

Analysing the data for the magnetic susceptibility in the high temperature phase we 
arrive at 

C+ = 0.1739(1) + 800 x (p c - 0.22165462) - 1.6 x (7 - 1.2372) - 0.0013 X (w - 0.832) (A14) 
and 

a Xi+ = -0.33(5)-150000x(/3 c -0.22165462) + 100x( 7 -1.2372)-1.3x(u;-0.832) . (A15) 

The corresponding result for the low temperature phase are 
C_ = 0.03695(2) - 200 x (/3 C - 0.22165462) - 0.35 x (7 - 1.2372) - 0.001 x (w - 0.832) (A16) 
and 

a Xi _ = -1.6(2) + 20000 x (fi c - 0.22165462) + 120 x (7 - 1.2372) - 7x(w- 0.832) . (A17) 



The ratio C+/C- = 4.706(8) is consistent with = 4.713(7) obtained in ref. |3jJ by 

analysing Monte Carlo data obtained for the Blume-Capel model at D = 0.655. Note that 
our estimates are slightly smaller than C + /C- = 4.78(3) obtained from series expansions 

FT 



3. The energy density 

In order to compute the thermodynamic Casimir force, we need the energy density of the 
bulk system for a large number of values. In princible one might use the results given in 



451 ]. However these estimates are not sufficiently accurate for our purpose. Instead, we have 
combined the analysis of the high [39[ and low [!o[ temperature series with the results of 
our Monte Carlo simulations discussed above. The analysis of the high temperature series 
is simpler and the results are more accurate than that of the low temperature one. This is 
due to the fact that the high temperature series converges up to the critical point, while this 
is not the case for the low temperature series. 

In the neighbourhood of the critical point the energy density behaves as 

Ebuik = E ns - C ns t + ... + ail^-^l + b±\t\ e + ...) + ... (A18) 

We analysed both series using differential approximants. In particular, we used the second 



order differential equation given in eq. (6.16) of ref. 46]: 

u 2 Q 2 {u)g"{u) + uQ^g'iu) + Q (u)g(u) = R{u) (A19) 

where Q2(u), Qi(u), Qo(u) and R(u) are polynomials in the expansion variable u of the 
order J, K, L and M, respectively. These polynomials are fixed by the requirement that the 
function g(u) has the correct expansion in u up to the highest known order. The differential 
eq. (IA19p is used, since it is known that its solution behaves as 

g(u) = g ns (u) + ai(«)(« c - u)~ Xl + a 2 (u)(u c - u)~ X2 (A20) 
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where g ns (u), a\{u) and 02 (u) are analytic functions. 

Usually one sets ^2(0) = 1. Therefore J + K + L + M = N — 2, where N is the order of 
the last known coefficient of the series. We biased the analysis by using our estimate (IA2j) 
of the inverse critical temperature and our estimates of v and u [7|. This way additional 
coefficients of the polynomials are fixed and one gets J + K + L + M = N + 3. For a 



detailed discussion we refer the reader to section 6 of ref. |46j. We solved the differential 
equation (|A19[) numerically by using the Runge-Kutta method. 

In the high temperature phase Arisue and Fujiwara [39J computed the free energy density 
of the bulk system as a series in v — tanh(/3) up to 0(v m ). Note that the coefficients 
of odd orders vanish and hence the free energy density can be expressed as a series in 



u = v Since we are aiming at the energy density, we have actually analysed 



E = . (A21) 

ou 



The energy density is then given by 



E bulk = -% = = 2tanh(/3)[l - tanh 2 (/3)] E (A22) 



d(3 du d(3 
The free energy density is given by 
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- f{fi) = In 2 + 3 ln(cosh(/3)) + ^ + 0(v 48 ) (A23) 

i=0 

where the coefficients ai are given in table I of the preprint version of ref. [3 

First we analysed E to obtain estimates of the energy density. We computed x 2 — 
Y^i[( E series(/3i)-E M c(/3i))/e(Pi)} 2 , where E series ((5i) and E MC ((5i) are the estimates obtained 
from the analysis of the series and from the Monte Carlo simulations, respectively, and e(/?j) 
is the statistical error of the Monte Carlo result at the inverse temperature We find that 
a large fraction of the possible choices of J, K, L and M result in a x 2 /d.o.f.^ 1.03. About 
91% of the possible choices have x 2 /d.o.f.< 1.073 and about 92.5% have x 2 /d.o.f.< 1.305. 

We computed numerically E ns , C ns , a + and a + b + as defined by eq. (1A18I) . Averaging 
over all choices of J, K, L and M with x 2 /d.o.f.< 1.073 we get 

E ns = 0.9906058(8) + 32 x - 0.22165462) 
- 0.0069 x (u - 0.63002) 

+ 0.0000072 x(u- 0.832) , (A24) 

C ns = -29.07(3) - 234000 x ((3 C - 0.22165462) 

- 1960 x (v- 0.63002) 

- 0.86 x (lo - 0.832) , (A25) 

a+ = -25.715(12) - 92500 x (p D - 0.22165462) 

- 1390 x (u - 0.63002) 

- 0.244 x (w - 0.832) (A26) 
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and 



a + b + = 3.87(28) - 1300000 x (ft c - 0.22165462) 
- 2900 x(v- 0.63002) 

+ 13 x (u - 0.832) . (A27) 

The number given in () is the variance over all choices of J, K, L and M with x 2 /d.o.f.< 
1.073. It might serve as a lower bound of the systematic error of the analysis of the series. 
Since the estimates for E ns and C ns obtained here are in good agreement with those obtained 
from the finite size analysis of Monte Carlo data given above, we are confident that also in 
the case of a + and a + b + the variance over the choices of J, K, L and M is a realistic estimate 
of the systematical error. Analysing the series for the free energy density itself we get 

- f ns = In 2 + 0.0847028611(4) + 0.99 x (ft c - 0.22165462) + 0.000001 x (u- 0.63002) (A28) 

The estimate of f ns strongly depends on the input value for ft c . The dependence on v is 
small and that on u can be ignored. 

In order to calculate the energy density that is needed as input to compute the thermody- 
namic Casimir force we have picked, to some extend ad hoc, the approximant characterised 
by J = 7, K = 7, L = 5 and M = 6 which is characterized by the fact that the order of 
all four polynomials is similar, x 2 /d.o.f= i- 029 and E ns = 0.9906063 for ft c = 0.22165462, 
v = 0.63002 and u = 0.832 fixed. Comparing with other acceptable choices for J, K, L and 
M we find that e.g. for ft = 0.2216 the differences are of the order 10~ 7 and for ft = 0.22 of 
the order 10~ 8 . Compared with the statistical error of [E(L + d/2, ft) — E(L Q — d/2, ft)]/d, 
see eq. ( 18"3"|) . errors of this size are negligible. 

In the low temperature phase, Vohwinkel [!o[ computed the energy density as a series in 
u = exp(— 4/3) up to 0(u 32 ). Unfortunately in this case there is no choice of J, K, L and M 
that allows to fit our Monte Carlo data down to ft = 0.2219. The best that we could find 
are the two choices J = 9,K = 6,L = 7 and M = 13 and J = 20, K = 6, L = 3 and M = 6 
that fit our Monte Carlo data with an acceptable x 2 /d.o.f. for ft > 0.228 and ft > 0.231, 
respectively. The linear combination 0.8155^9^,7,13 + 0. 1845-^20,6,3,6 fits all of our data in 
the low temperature phase with x 2 /d.o.f.= 1.25. 

Since this result is not fully satisfying, we fitted our data with various ansaetze based on 
eq. (I A 1 8 [) . In particular the ansatz 



E = E ns — C ns t + d ns t 2 + a4-ty- a + a_6_(-t) 1 -° +e + b(-tf- a + c(-tf- a - d (A29) 

fits our data up to ft = 0.246 with x 2 /d.o.f.= 1.15, where we have fixed E ns = 0.9906065, 
C ns = -29.07, a = 0.10994 and u = 0.832. Fitting all 55 data points up to ft = 0.241 we 
get for the free parameters a_ = 47.9436, a_6_ = —16.336, b = —363.5, d ns = 269.2 and 
c = 287.3. In order to calculate the bulk energy that is needed for the computation of the 
thermodynamic Casimir force we used for ft > 0.228 the linear combination 0.8155-^9^,7,13 + 
0.1845£ , 2 o,6,3,6 °f approximants and for 0.228 > ft > ft c we used eq. (IA29|) together with the 
results for the free parameters quoted above. For a quite large range of ft the two approaches 
to represent the bulk energy give consistent results. For 0.2219 < ft < 0.2394 the difference 
between the two is less than 3 x 10~ 6 . The deviation of our result from that of ref. j45| is 
typically of the order 10 -5 . 
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Taking into account various fits and in particular computing the dependence of the result 
on the values of the input parameters, we arrive at 

a_ = 47.96(1) + 2350000 x (/3 C - 0.22165462) 
+ 2500 x (u - 0.63002) 
-0.16 x (a; - 0.832) 

- 0.44 x {C ns - 29.1) 

- 3700 x (E ns - 0.9906065) (A30) 

and hence 

4 1 = -— = 0.5362(20) (A31) 

A_ a_ 

which is fully consistent with the estimate A + /A_ = 0.536(2) obtained by studying the 



Blume-Capel model at D = 0.655 [311 ] . Note that the error of our estimate of A + /A_ is 
dominated by the uncertainty of C ns that we use as input for our fits in the low temperature 
phase. Here we have taken the error of the estimate obtained from the finite size scaling 
analysis at the critical point, eq. ( 1A6I) . The systematic error of the estimate obtained from 
the analysis of the high temperature series is likely smaller, but difficult to estimate. The 



authors of [44| quote A + /A_ = 0.530(3) which is slightly smaller than our results. For a 



summary of estimate presented in the literature see table IV or ref. 44 . 
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